Counting attractors in synchronously updated random Boolean networks 
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Despite their apparent simplicity, random Boolean networks display a rich variety of dynamical 
behaviors. Much work has been focused on the properties and abundance of attractors. We here 
derive an expression for the number of attractors in the special case of one input per node. Ap- 
proximating some other non-chaotic networks to be of this class, we apply the analytic results to 
them. For this approximation, we observe a strikingly good agreement on the numbers of attractors 
of various lengths. Furthermore, we find that for long cycle lengths, there are some properties that 
seem strange in comparison to real dynamical systems. However, those properties can be interesting 
from the viewpoint of constraint satisfaction problems. 

PACS numbers: 89.75.Hc, 02.10.Ox 



INTRODUCTION 

Random Boolean networks have long enjoyed the at- 
tention of researchers, both in their own right and as 
simplistic models, in particular for gene regulatory net- 
works. The properties of these networks have been stud- 
ied for a variety of network architectures, distributions 
of Boolean rules, and even for different updating strate- 
gics. The simplest and most commonly used strategy is 
to synchronously update all nodes. Networks of this kind 
have been investigated extensively, see e.g. 

The networks we consider arc, generally speaking, such 
where the inputs to each node are chosen randomly 
with equal probability among all nodes, and where the 
Boolean rules of the nodes are picked randomly and inde- 
pendently from some distribution. The network state is 
updated synchronously, with the state of a node at any 
time being a function of the state of its inputs at the 
previous time step. 

In this work we determine analytically the numbers 
of attractors of different lengths in networks with con- 
nectivity (in-degree) one. We compare these results to 
networks of higher connectivity and find a remarkable de- 
gree of agreement, meaning that networks of single-input 
nodes can be employed to approximate more complicated 
networks, even for small systems. For large networks, a 
reasonable level of correspondence is expected. See 
on effective connectivity for critical networks, and on 
the limiting numbers of cycles in subcritical networks. 

Random Boolean networks with connectivity one have 
been investigated analytically in earlier work 0. 
In those papers, a graph-theoretical approach was em- 
ployed. In this work, we use another approach |lfll | which 
is based on the ideas of fixed point analysis. Our ap- 
proach is a powerful tool for counting the average num- 
bers of attractors, but it does not give any direct infor- 
mation on the properties of the attractor basins. 

For long cycles, especially in large networks, there are 
some artefacts that make comparisons to real networks 



difficult. For example, the integer divisibility of the cy- 
cle length is important, see e.g. Also, 
the total number of attractors grows superpolynomially 
with system size in critical networks [Sl llOj , and most of 
the attractors have tiny attractor basins as compared to 
the full state space. In this work, such artefacts become 
particularly apparent, and we think that long cycles are 
hard to connect to real dynamical systems. 

Comparisons to real dynamical systems, on the other 
hand, still seem to be relevant with regard to fixed points 
and some stability properties 0, Il2l | . An interesting way 
to make more realistic comparisons regarding cycles is to 
consider those attractors that are stable with respect to 
repeated infinitesimal changes in the timing of updating 
events 0. 

Furthermore, the problem of finding attractors with 
small attractor basins, for a given network, can be seen 
as a hard constraint satisfaction problem. Those prop- 
erties that are artefacts when comparing to dynamical 
systems, could be keys to the understanding of real life 
combinatorial optimization problems. 

Throughout this paper, N denotes the number of nodes 
in the network, and L the length of an attractor, be it a 
cycle (L > 1) or a fixed point (L = 1). For brevity we 
use the term L-cycle, and understand this to mean an 
attractor such that taking L time steps forward produces 
the initial state. When L is the smallest positive integer 
fulfilling this, we speak of a proper L-cycle. We denote 
the number of proper L-cycles with Cl- The average 
over networks of a certain size is {Cl)n, so the average 
number of network states that are part of a proper L- 
cycle is L{Cl)n- Related to this is {^l)n, the number 
of states that reappear after L time steps and hence are 
part of any L-cycle, proper or not. 

To calculate the average number of attractors, we use 
the same approach as in [loj . whereby we first trans- 
form the problem of finding an L-cycle into a fixed point 
problem, and then find a mathematical expression for the 
average number of solutions to that problem. 
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In the case that every node has one input, {flL)N 
can be calculated analytically for any A''. The limit 
{^l)oc has a relatively simple form, and this limit can be 
mapped to the 1-input case for any subcritical network of 
multi-input nodes, as shown in Furthermore, critical 
networks of multi-input nodes seem to show strong sim- 
ilarities to the corresponding networks of 1-input nodes. 



THEORY 

Assume that a Boolean network performs a proper L- 
cycle. Then, each node performs one of m = 2^ series 
of output values. We call these L-cycle series. We as- 
sociate each L-cycle series with an index i, such that 
i E {0, 1, . . . , TO — 1}. For convenience, we let the con- 
stant X-cycle series have the indicies and 1, in such a 
way that a constant false output has index whereas a 
constant true output has index 1. 

If we view each L-cycle series as a state, an L-cyclc of 
the entire network turns into a fixed point (in this en- 
larged state space). L{Cl)n is then the average number 
of input states (choices of L-cycle series), for the whole 
network, such that the output is the same as the input. 

It is inconvenient to work directly with {Cl)n- Let 
{^II) denote the average number of states that reappear 
after L timesteps. Such a state is part of a proper i- 
cycle where £ is a divisor to L. {Cl)n can be calculated 
from (f2i)jv, using the set theoretic inclusion-exclusion 
principle. 



Expressing {C'l)n in terms of {^Il)n 

Let Cl denote the set of network state sequences 
that represent proper L-cycles. Similarly, let ujl de- 
note the non-proper counterpart to Cl, meaning that 
ujl = [J^^Ce where £\L means that £ divides L. 

Consider a given network, and let Q denote the set 
of sequences of states, Q, that are consistent with the 
network dynamics. Let i2>l = tOL Q and Cl ~ Cl Q- 
Then, Ql = I'^lI, and the number of proper L-cycles of 
that network is given by Cl = ^ClI- To calculate Cl, 
we start by expressing Cl in terms of cue. We sec that 



Cl^u;l\ U 



UJi = UJl 



l<i<L 
l\L 



\ U ^L/d^ 
d prime 
d\L 



(1) 



because any positive £ dividing L is also a divisor to a 
number of the form L/d, where d is a prime. 

Then, the inclusion-exclusion principle, applied to 
eq. 1^, yields 



se{o,i}'^L 



(2) 



where s = j:ti ^^(s) = UtiidlY' and 4, 
are the prime divisors to L. For averages over randomly 
chosen A'^-node networks, we get 



se{o,i}'^L 



-^y{^L/dL{s))N 



(3) 



Basic expressions for {Q,l)n 

Consider what a rule does when it is subjected to L- 
cycle series on its inputs. It performs some Boolean oper- 
ation, but it also delays the output, giving a one step dif- 
ference in phase for the output L-cycle series. Let ^^(x) 
denote the probability that a randomly selected rule will 
output L-cycle series i, given that the input series are se- 
lected from the distribution x = (xq, . . . ,Xm-i)- Then, 
we get 



neN'" 

n=N 



m—1 



(4) 



where n = uq + ■ ■ ■ + Um-i- This is the same expression 
as presented in 0. 

Let eo be the vector with elements Cq = 5io, where 6 
is the Kronecker delta. For networks of 1-input nodes, 
A^(x) is affine and can be expressed as 

Al(x) = (x-eo)-VAl+^l(eo) (5) 
= X • VA^ + dioco + Siici , (6) 

where r'~^ and arc the fractions of nodes that copy and 
invert their input, respectively, while cq and ci arc the 
fractions of constant nodes that output false and true, 
respectively. The elements djA'^ of the gradient WA]^ arc 
nonzero only if the L-cycle series i can be the output of 
a node receiving the L-cycle series j on its input. This 
is true if the series j is the series i, or the inverse of i, 
rotated one step backwards in time. Let (t>^{i) and (/'^(i), 
respectively, denote those values of j. Then, 

Let Bl be an TO X TO matrix with elements 

B'l = djA^ + 5,o(co - 1) + ci5,i ■ (8) 

Then, 

AL(x)=BLX + eo (9) 
for all X such that xq + . . . + x„i-i = 1. Furthermore, 



^(BlX + Gq), = 1 



(10) 



1=0 



for all X. This property will be important later on. 
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Prom the theory of Hnear algebra, we know that any 
matrix can be transformed to a triangular matrix. Thus, 
B/, can be rewritten as 



(11) 



where Cl is an invcrtable matrix and is a triangular 
matrix. As we shall see, it turns out to be sufficient to 
know det(l — CDl) as a function of C to calculate (^Jl) n- 
Furthermore, 

det(l - CDl) = dct(l - CBl) (12) 

because 

(1-CBl) = C^1(1-CDl)Cl . (13) 
Thus, we want to calculate dct(l — CBl). 

The structure of Bl 

To understand the structure of B^, we consider sub- 
spaces spanned by i-cycle series indices of the type 

{«, <^L (i) : (0 > ° (*). </'l ° "^L («)'■■•} > containing aU 
possible results of repeatedly applying (p1 and to i. 
We call those sets invariant sets of L- cycle series, which 
is the same as the invariant sets of L-cycle patterns in 
[lo| , but formulated with respect to i-cycle series instead 
of L-cycle patterns. Let p^,! • ■ • i Pl^~^ denote the invari- 
ant sets of i-cycle series, where H]^ is the number of such 
sets. For convenience, let p° be the invariant set |0, Ij. 
These definitions are consistent with the notation in [lOj , 
meaning that is the same number and po corresponds 
to the same invariant set. 

The elements of B^ are given by 

= + + S^oico - 1) + Saci , (14) 

which means that B^ and 1 — CBl are block triangu- 
lar matrices, with blocks corresponding to the invariant 
sets. Consequently, det(l — CBl) is the product of the 
determinants of all blocks on the diagonal. 

Let r = r'~^ +r^ and Ar — r'~^ — r^. Then, the determi- 
nant of the block corresponding to is 

V -C{r^ + ci) 1-C(r^+Ci)y' 

(15) 

because r'^ + r^ + cq + ci ~ 1 . 

To calculate the determinants of the blocks corre- 
sponding to the other invariant sets, we need to explore 
the structure of the invariant sets. Consider an invariant 
set of L-cycle series, p^. Let £ be the length of p'l, mean- 
ing that i is the lowest number such that, for i € p'l, 
{4>2Y{'i) is either i or the index of series i inverted. If 
{(j)^Y{i) ~ i, wc sa.y that the parity of p'l is positive. Oth- 
erwise the parity is negative. The structure of an invari- 
ant set of L-cycle series is fully determined by its length 



and its parity. Such a set can be enumerated on the form 
{0g(O,...,(0g)^(*),0i(z),...,</.i 



(z), } and for 
b2Y-\i) = t ior 



positive parity (</>/;) («) = i, while 0^ o | 
negative parity. 

Let strings of T and F denote specific L-cycle series. 
Then (j>'^{FFFT) = fftf and (/)^(ffft) = ttft. Exam- 
ples of invariant sets of 4-cycle series are {ffft, fftf, 

FTFF, TFFF, TTTF, TTFT, TFTT, FTTT} and {fTFT, 

tftf}. The first example has length 4 and positive par- 
ity, while the second has length 1 and negative parity. 

Let and denote the blocks in the diagonal of 
B L that correspond to invariant sets of length £ with pos- 
itive and negative parity, respectively. Before we present 
the form of R^ for general £, we consider a few examples. 

The invariant set {ftft, tftf} corresponds to 



Rr 



(16) 



where the first row corresponds to the 4-cycle series ftft 
while the second row corresponds to tftf. Correspond- 
ingly, we define 



R^ 



+ — 



(17) 



Using this definition, the other example of an invari- 
ant set, {ffft, fftf, FTFF, TFFF, TTTF, ttft, TFTT, 

fttt}, corresponds to 



Rt 



/ 




VR+ 



R+ 








R+ 







R+ 

/ 



(18) 



with rows and columns connected to the 4-cycle series 

FFFT, TTTF, FFTF, TTFT, FTFF, TFTT, TFFF, and FTTT, 

in that order. 

For general £ and parity, we get 



/ R+ 




R» 



VRf 










R+ 








\ 




R+ 
/ 



(19) 



where order of the rows (and columns) is given by in- 
dex sequences of the type (j)^{i), 4'\{i), ■ • ■, {<l)^Y{i), 4^ ° 

The unitary tranformation matrix 
Si 

transforms Rj*^ according to 

Ri ^SiR^Si- 



1 

71 



1 -1 
1 1 



(20) 



(21) 
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Let 



/Si 






Si 



\ 



0\ 





Si/ 



(22) 



Then, 
det(l - CR-?) 



det[S/(l-CR£^)Sf] 
det(l 



det 



1 




-CR+ 
1 








= [l-(rCn[lT(ArCn 





-CR+ 






\ 


<R+ 

1 / 



(23) 
(24) 



(25) 



(26) 



To calculate det(l — CBl), we need to find the dis- 
tribution of lengths and parities of the invariant sets of 
i-cycle series. If an L-cycle series belongs to an invariant 
set with length £ it will be itself or itself inverted after £ 
timesteps, depending on whether the invariant set is of 
positive or negative parity. This gives a periodicity of £ 
or 2£ timesteps. Hence, for L-cycles, there will be invari- 
ant sets of length £ and positive parity if £\L, whereas 
sets of negative parity are present if 2£\L. 

If invariant sets of a specific length and parity are 
present, their number is independent of L, because the 
basic form of the series in such sets does not change with 
L. Only the number of basic repetitions differs. Let 
and denote the numbers of invariant sets of length £ 
with positive or negative parity, respectively. Then, 



det(l - B^C) = n [1 - [1 - 



xn[i-K)r^ [i+(A<)T^ > 

(27) 

where the Kronecker delta, 5n , takes care of the special 
case of the block of , corresponding to pj, • This block 
has the determinant 1 — ArC, instead of the (1 — rC,){l — 
ArC) it would have if it obeyed eq. H26|l . 



Calculation of {Ql}n via tensor calculus 

Here, we derive how (r2i)jv can be determined from 
dct(l — BiC)- From eq. 10} and eq. jnj, we see that 



nGN'" 



(28) 



Eq. (|28|l can be rewritten, using (m -|- m)-dimensional 
tensors such that each element is indexed by two of the 
vectors k. A, /2, e N™. Each such vector corresponds 
to a distribution of L-cycle patterns, and k. A, fi and v 
denote the sums of the elements in each vector. We get 



LIN 



Tr(Gi) , 



where 



{GLYa = Sf^N 



N 



llliBLiy/N + eoY 



3=0 



and the trace operator is defined as 

Tr(X) ^ E 4 ' 

0</jn....,/^m_i 

for an (m -I- m)-dimensional tensor X. 

For convenience, we define the operation n— as 

k~l 



(29) 



(30) 



(31) 



(32) 



We also choose the convention to interpret nP as an empty 
product, independently of the value of n. This means 
that O'' = 1, and division must be treated with special 
care. Let 



and 



m — 1 

- n 

3=0 



3=0 



(33) 



(34) 



M is triangular in the sense that kq < /ip, . . . , k,„_i < 
/i„i-i for all non-zero elements Mj^. Hence, M has an 
inverse that obeys fio < kq, . . . , fim-i < k„i-i for all 
non-zero elements (M~^)^. 
Letting M act on Gl yields 



m— 1 

X l[[{BLi^/N + eo), 

3=0 

771—1 



N -K 



(35) 



3=0 
ni—1 

3=0 

rn—1 

3=0 



(36) 
(37) 
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because X^j^o i^L^ + ^o)j = 1- Eq. can be rewrit- 
ten to 

mGLfp = T7^ n [(Bl^). + NS,or (38) 



i=o 



m— 1 



(39) 



Let Z5(X), for an arbitrary m x m matrix X, denote 

the coefficients of the formal expansion of IIJl^^ [(^2)j]'^^ 
in such a way that 

m— 1 m — 1 

n[(xz),f^ ^ ^ zf(x) n ■ (40) 

j=0 0<Ko,...,fOm-l i=0 

Then, we see that the operator Z satisfies the condition 

Z(XY) = Z(X)Z(Y) (41) 

for all m X m matrices X and Z. Furthermore, Z(X) is 
block diagonal for any X, in the sense that k = X for all 
nonzero elements Z£(X). 
Define djjp according to 



J=0 



(42) 



for arbitrary m-vectors v and /i. Now, eq. H39() can be 
rewritten as 



where 



and 



MGl = NEZ(Bl)M 



NS = S 



(43) 



(44) 



(45) 



The block diagonal structure of Z(X) yields 

NZ(X) = Z(X)N (46) 

for all X, because N is diagonal with diagonal elements 
that are constant for constant k. According to eq. 
B L can be written as 



where is triangular. Thus, 

MGl = NEZ(C^1DlCl)M (48) 
Z(CL)MGL = NELZ(Di)Z(Ci)M (49) 
Z(Ci)MGiM-iZ(C^i) = 

= NElZ(Dl)Z(Cl)MM-1Z(C£^) (50) 

where El = Z(Cl)EZ(C^^). Because M^^ is triangu- 
lar with the lower indicies (acting to the left) less than 
or equal to the upper indicies, right multiplication with 
always yields a convergent result. Similarly, Z{Cl) 
and Z(C^^) are well-behaved with respect to both left 
and right multiplication, as a consequence of their block 
diagonal structure. Thus, both sides of the equality in 
eq. 1501 are well-defined. 

Let Tl = Z(Cl)MM-iZ(C£^), which yields 

Z(Cl)MGlM-1Z(C^1) = NELZ(Di)Ti . (51) 

The tensor MM"^ tells how to express moments in terms 
of combinatorial moments. Each moment can be ex- 
pressed as a sum of the combinatorial moment of the 
same order and a linear combination of combinatorial 
moments of lower order. Hence, 



Me(Af-M;: = forA>K. 

K. V / f-L kX — 



(52) 



This property is conserved as MM ^ is transformed 
by Z{Cl), because 

Z§{Cl) = for AC ^ , (53) 

meaning that 

= % for fi>iy. (54) 

This is also true for E^, i.e., 

eI = eI = d.~^ for A > K . (55) 

Because T^, and E l are triangular (with the same ori- 
entation) with unitary block diagonal, and N and Z(Di) 
are block diagonal, we get 

{nL)N = Tr(GL) = Ty[NElZ{T>l)Tl] = Tr[NZ(DL)] . 

(56) 

The triangular structure of yields that the diagonal 
elements of Z(Di) are given by 



?n — 1 

z5(Di)=n(^f)'^^ 

J=0 



(57) 



Thus, 



Tr[NZ(D,)] = J^U^:^ (58) 

0<K.o,...,Km-l j — 



m—1 



(47) 
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and 



N 



N- 



m — 1 



LIN 



u—{) 0<ko,...,h:tti-i J— 

TV 



(59) 



E 

v=0 



N'-^ 1 (F 



771 — 1 

E lii^LO"^ (60) 
C=0 0<Ko, j=o 



Km-l 

ni~ 1 



1 d 



TV 



n — 1 



(61) 
(62) 



Because 



n (1 - D^iO = dct(l - DiC) 
= det(l - BiC) 



we get 



(f2L)Ar = 1 + 



^_d_ 

NdC 



N 



det(l - BlC) 



(63) 
(64) 

(65) 



where 



det(l - B^C) = n [1 - " [1 - 

2l\L 



(66) 



Calculation of 



Let ^/i^ and ■0^, respectively, denote the sets of (infi- 
nite) time series of true and false, such that each series 
is identical to, or the inverse of, itself after i time steps. 
Then, the set of time series that are part of an invariant 
set of L-cycles with length i and negative parity, , is 
given by 



Jf ^YY\ U 



(67) 



d odd prime 

d\£ 



For positive parity, we get 



prime 
d\£ 



The numbers of elements in are given by 2iJ^, 
where are the numbers of invariant sets with length 

(69) 



Then, the inclusion-exclusion principle yields 
1 

se{o,i}''''' 



where s = J2lLi d{s) = Hl^lil^f)''' ^^'^ d], . . . , d^ are 
the odd prime divisors to 1. Similarly 

Ji-hT. E (-l)^°+^2^/''^(^-^) - i J,-/^ , (70) 
soe{o,i} se{o,i}*'* 

where d(so:S) = 2''°d(s) and J^"^, = if 1/2 is not an 
integer. Insertion of Eq. into Eg. 1701 vields 

Jt^ TiY. E (1 + so)(-l)^»+^2^/'^^(^»'^) , (71) 
soe{o,i} se{G,i}^i 



which also can be written as 



(72) 



Calculation of (p.h)N for small and large A'^ 

Iterative application of cq. (|72|) and the identity 

[1 - (ArC)^] [1 + (ArC)^] = 1 - (ArC)^^ (73) 
to cq. H66|) yields 

det(l - BiC) = 



n [l-(rC)T' -'^/=-*"[l-(ArC)^]' 



Ljt odd 



n 



(74) 



Lji even 



where J ^ is given by cq. iPll . 

The simplest way to calculate {^l)n is to work with 
power series expansions. If 



{nL)N is given by 



oo ^ 



(l-BiC) ' 



(75) 



k=0 



(76) 



where ,2 is the empty set if i/2 is not an integer. 



where the operation A^- is defined in cq. (|32|) . 
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For r < 1, C = 1 is within the convergence radius of 
the power scries in cq. (I75|) . Hence, we get 



hm 1 + —-7- 

N->oo\ N d( 



N 



1 



1 



det(l - BlO dct(l - Bl) ' 

(77) 

which means that the large N hmit of {^II) n is given by 



1 



L/oo 



det(l - Bi) 



(78) 



This resuh is consistent with the calculations in Q , which 
also are valid for subcritical networks with many inputs 
per node. 

For r = 1, the dominant contribution to {^l)n comes 
from the pole at ^ = 1. For |Ar| < 1 and large iV, we 
get 



N 



IL 



c=o(i-C) 



(79) 



where is the total number of invariant L-cycle sets 
and 7l is a constant. For large N ^ a power series expan- 
sion of - C)^"-"^] yields 



N 



k=0 
/•oo 



fHL -2 + k 
TV^V k 



IL 



(80) 
(81) 



where 7^ is a constant. The dominant terms in eq. H8U|) 
satisfy H]^ <^ k ^ N in the large N limit, and Stirling's 
formula has been applied to those. From cq. (|81|l . we see 
that the asymptotic behavior of (O^) ^ is the scaling 



(82) 



for critical r ~ 1 networks with |Ar| < 1. This also 
means that 



(83) 



for large N . 

If |Ar| = 1. the network only consists of copy oper- 
ators (if Ar = 1) or inverters (if Ar = — 1). Then, Hl 
should be replaced by the number of sets of L-cycle series 
that are invariant under the present choice of operators. 
This number will be as least as large as Hl, because the 
invariant set can split, but not merge, when one operator 
is removed. 

One can also use complex analysis to retrieve the op- 
erator identity 



N 



(84) 



where e is a positive constant small enough to keep all 
poles of the given function outside |C| < £• This is es- 
pecially useful to calculate the total number of states in 
attractors; (r2o)7v- 

Eq. ()74|l can be rewritten as the exponential of a power 
series expansion. The expansion 



Mi-^) = E 



(85) 



yields 

-lndet(l - BiC) 



l\L 
Ljl odd 



2t 



l\L 
L/^ even 



(86) 



On- 



where J^°^^ = J^-J'j^-Sn and J^°™" = 2.J^-J- _ 

Let (a, b) denote the greatest common divisor of a and 
h. Then, the substitution k = j£ and reordered summa- 
tion gives 

-lndet(l-BLC) = 



fe=l e\{k,L) 
Ljl odd 



(87) 



E 



f|(fe,L) 
hli even 



E [-r'= + (Ar)^]J-^ 



fc>0 i\{k.L) 
L/{k,L) odd {k,L)/e odd 

+ E^E^r 



fc=i 



l\{k,L) 



The inner sums can be calulated using the connection 
between J 7 and ip7 . Thus, 



Ak 

Indet(l-BiC) = E ^{-r'' + iM'']'^ 



-)(k,L)-l 



k>0 
L/{k,L) odd 



E 

k=l 



(2 



(fe,L) 



(89) 



For L = 0, we get 

00 

-lndet(l -BoC) = E 



k=l 



(2rC)^ (rC) 



fe-i 



(90) 
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which means that 



L-cycles 



1 _ 1-rC 

det(l - BoC) ~ 1 - 2rC 



and 



N 



dc 



^NQ 1 - rC 



(91) 



(92) 



The A^-dependent part of the integrand in eq. (|92|l has 
a saddle point atC = l. Ifr<l/2, the path of steepest 
descent integration does not include any other poles than 
the desired one at C = 0. On the other hand, if r > 
1/2 the pole at C = l/(2r) gives another contribution to 
(r2o)Af- The borderline case, r = 1/2, is similar to the 
case r = 1 for {^l) n with L 7^ 0. 

Considering the properties mentioned above, we find 
that the asymptotic behavior of (fio) n is given by 



1 - r 
1 - 2r 



for r < 1/2 

{%)n ~ { hVW r = 1/2 . (93) 

pre(ln2r-l + l/2r)W for r > 1/2 

RESULTS 

Our main results from the analytical calculations are 
the expressions that yield (^lL)N and {Cl)n for finite 
TV, and their asympotic growth. See eqs. ©, (fT^Tzfljl 
and l|91() on expressions for general and eqs. (|78|1 . 
(IS^ . and (|^ on expressions valid for the high N 
limit. Fig. n shows the numbers of attractors with short 
length as a function of system size, for different but 
with Ar = 0. For critical networks, with r = 1, the 
asymptotic growth of (Cl) at is a power law, while (Cl) at 
approaches a constant for subcritical networks as N goes 
to infinity. 

For networks with Ar =^ 0, the prevalences of copy 
operators and inverters are not the same. Cycles of even 
length are in general more prevalent then cycles of odd 
length. An increased number of inverters strengthens 
this difference, while a low fraction of inverters makes the 
difference less pronounced. See fig. 12 which shows the 
symmetric case Ar = and the extreme cases Ar = ±r. 

The total number of attractors, {C)n, and the total 
number of states in attractors, (r2o)Af, can diverge for 
large N, even though the number of attractors of any 
fixed length converges. This is true for subcritical net- 
works with r > 1/2. See figs. 13 and The growth 
of {^Iq) N is exponential if r > 1/2. Interestingly, there 
is no qualitative difference in the growth of (fio) n when 
comparing the critical case, r = 1, to the subcritical ones 
with 1 > r > 1/2. 
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FIG. 1: The average number of proper L-cycles as a function 
of A*' for different L, for networks with single-input nodes, 
r = 1 in (a), and r = 3/4 (solid lines) and r = 1/2 (dotted 
lines) in (b). In (a), L is indicated in the plot. In (b), L is 
3, 5, 7, 1, 2, 4, 6, 8 for r = 3/4 and 7, 5, 3, 8, 6, 4, 2, 1 for 
r — 1/2, in that order, from bottom to top along the right 
boundary of the plot area. In (b), the curves for L = 3 and 
L = 5 for r = 3/4 essentially coincide at the right side of 
the plot, whereas they split up to the left, with L = 3 as the 
upper curve there. 



For r < 1/2, both {C)n and {ilo)N converge to con- 
stants for large N. In the border line case r = 1/2, 
(fio) AT diverges like a square root of N, but (C) n seems 
to approach a constant. See fig. QJ). 

The quantity (r2o)w has an interesting graph theoret- 
ical interpretation. Let A^activc be the fraction of nodes 
that are part of a loop of non-constant nodes. Then, 



(fio)jv = (exp(7Vactivoln2)) 



(94) 



which means that ln(exp(iVactive)) grows linearly with 
if 1/2 < r < 1. This stands in sharp contrast to 
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FIG. 2: The average number of proper L-cycles for networks 
with N = 100 and r = 3/4, as function of L. Ar = 3/4 (thin 
solid line), Ar = (thick solid line) and Ar = —3/4 (dotted 
line). Note the importance of what numbers divide L. 
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FIG. 3: {Qo)n (solid lines) and {C)n (dotted hues) for r = 
1/2, r = 3/4 and r = 1, in that order, from the bottom to 
the top of the plot. Ar = in all cases. Note that {^Io}n 
is independent of Ar. Due to the double logarithmic, {Cl)n 
appears to approach (Ql)n for r — 3/4 and r — 1. This is 
not the case — have in mind that a constant translation in a 
double logarithmic scale corresponds to a rise to a constant 
power. 

(-^active), which grows like \/N for r = 1 and approaches 
a constant for r < 1 as ^ oo , sec . This means that 
the distibution of A^active has a broad tail if r > 1/2. 

Like (r2o)jv is connected to the the number of nodes 
in active loops, {Ci)n for Ar = r is connected to the 
number of active loops. In this case, however, there is no 
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FIG. 4: The number of attractors with lengths L < Lmax 
in networks with A'^ single-input nodes, for different values of 
N. In (a) N = 10, 10^ . . . , lO'' for r = 1 (thin solid lines) 
and AT = 10, ... , 10* for r = 3/4 (thin dotted lines). In (b) 
iV = 10, 10^ 10^ (thin solid lines) for r = 1/2 and iV = 10 for 
r = 1/4 (thin dotted line). For all cases, Ar — 0. The thick 
lines in (a) and (b) show the limiting number of attractors 
when N — » oo. The arrowhead in (b) marks this limit for 
10^ for r = 1/2. The small increase in the number 
of attractors when Lmax is changed from 10'^ to 10^ indicates 
that (C)jv converges when A*' oo. Note the drastic change 
in the j/-scale between the case r > 1/2 and r < 1/2. 



striking difference between the average number of fixed 
points and the same quantity for a typical network. 

All the properties above are derived and calculated for 
networks with one input per node, but they seem to a 
large extent to be valid for networks with multi-input 
nodes. For such networks, we define r and Ar as mea- 
sures of disturbance propagation according to the follow- 
ing procedure: 

Find the mean field equilibrium fraction of nodes that 
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FIG. 5: Comparison between simulations for power law in- 
degree networks of size A'' = 20 (bold lines) and the corre- 
sponding networks with single-input nodes (thin lines). The 
fitted networks have identical values for r, Ar, and A'^. The 
solid lines show the number of fixed points, whereas the 
dashed and dotted lines show the numbers of 2-cycles plus 
fixed points and the total numbers of attractors, respectively. 
The probability distribution of in-degrees satisfies pK oc K 
where K is the number of inputs. The power law networks 
use the nested canalyzing rule distribution presented in Q. 

have the value true. Pick a random state from this 
equihbrium as an initial configuration. Let the system 
evolve one time step, with and without first toggling the 
value of a randomly selected node. The average fraction 
of nodes that in both cases copy or invert the state of 
the selected node are r'~ and r^, respectively. Finally, let 
r ~ r'~^ + and Ar = r'~^ — r^. For a more detailed 
description, see 0. 

From 0, we know that for subcritical networks the 
limit of {Cl)n as N — > oo is only dependent on r and 
Ar. Hence, we can expect that {Cl)n for a subcriti- 
cal network with multi-input nodes can be approximated 
with {Cl)'j^, calculated for a network with single-input 
nodes, but with the same r and Ar. For the networks in 
, with a power law in-degree distribution, this approx- 
imation fits surprisingly well. See fig. [S] 

For the critical Kauffman model with in-degree 2, one 
can do a similar comparison. The number of nodes that 
are non-constant grows like N'^^^ for large N. See |3.llo|. 
Furthermore, the effective connectivity between the non- 
constant nodes approaches 1 for large N j^. Hence, one 
can expect that this type of iV-node Kauffman networks 
can be mimicked by networks with N' = iV^^^ 1-input 
nodes. For those networks, we choose r = 1 and Ar = 0, 
which are the same values as for the Kauffman networks. 

For large N, {Cl)n in the Kauffman networks grows 
like n'--^^~^^^^ , where is the number of invariant L- 
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FIG. 6: Comparison between critical K — 2 Kauffman net- 
works (thick lines) and the corresponding networks of single- 
input nodes (thin lines). The size of the single- input networks 
is set to N' = N^^'^. r = 1 and Ar = 0, consistent with 
the Kauffman model, (a) The number of proper L-cycles for 
the L indicated in the plot. For the Kauffman networks, the 
numbers have been calculated from Monte Carlo summation 
for those network sizes where could could not be calculated 
exactly. See p^ . The number of fixed points is 1, indepen- 
dently of A'^, for both network types, (b) Total number of 
attractors. This quantity has been estimated by simulations 
for the Kauffman networks, using 10^ 10^ 10"* and 10^ ran- 
dom starting configurations per network. 



cycle sets. For the selected networks with 1-input nodes, 
we have {ClYn « iV'^^^-i)/^ cx iV^^^-i)/^ for large N' , 
see eq. |HS1 This confirms that the choice N' = N^^-^ is 
reasonable, but it does not indicate whether the propor- 
tionality factor in N' oc N^/^ is close to 1. This factor 
could also be dependent on L, as can be seen from the 
calculations in [lOj. However, this initial guess turns out 
to be surprisingly good, see fig. 05,. 
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Prom the good agreement for short cycles, one can ex- 
pect a similar agreement in the total number of attrac- 
tors. This is investigated in fig. |H1 For networks with 
up to about 100 nodes, the agreement is good, and the 
extremely fast growth of (Cl)^ for larger N is consistent 
with the slow convergence in the simulations. 

SUMMARY AND DISCUSSION 

Using analytical tools, we have investigated random 
Boolean networks with single- input nodes. We extract 
the exact distributions of cycles with lengths up to 1000 
in networks with up to 10^ nodes. Furthermore, we find 
some interesting scaling properties that hold for large N. 

As has been pointed out in earlier work we see 
that a small fraction of the networks have many more 
cycles than a typical network. This property becomes 
more pronounced as the system size grows, and changes 
the scaling of the average number of states belonging 
to cycles drastically. For networks that have the stabil- 
ity parameter r > 1/2, the average number of states in 
attractors grows exponentially with the system size 
whereas this number grows exponentially with ^/N for 
a typical network if r = 1, and approaches a constant if 
r < 1. 

The dynamics in random Boolean networks with multi- 
input nodes can to a large extent be understood in terms 
of the simpler single-input case. In a direct comparison 
to critical Kauffman networks of connectivity two and to 
subcritical networks with power law in-degree, the agree- 
ment is surprisingly good. 

Our results highlight some previously observed arte- 
facts in random Boolean networks. The synchronous up- 
dates lead to dynamics that largely is governed by integer 
divisibility effects. Furthermore, when counting attrac- 
tors in large networks, most of them are found in highly 
atypical networks and have attractor basins that are ex- 
tremely small compared to the full state space. 

In hj, a new concept of stability in attractors of 
Boolean networks is presented. To only consider that 
type of stable attractors is one way to make more relevant 
comparisons to real systems. Another way is to focus on 
fixed points and stability properties as in [^Il2l|. Further- 
more, the limit of large systems may not always make 
sense in comparison with real systems. Small Boolean 
networks may tell more about real systems than large 



networks would. 

Although there are problems in making direct compar- 
isons between random Boolean networks and real sys- 
tems, we think that insight in the dynamics of Boolean 
networks will improve the general understanding of com- 
plex systems. For example, can real systems have lots 
of attractors that are never found due to small attractor 
basins, and what implications would such attractors have 
on the system? 

Another interesting viewpoint is to compare with com- 
binatorial optimization problems, e.g. scheduling and 
digital circuit design. The question of finding an at- 
tractor in a random Boolean network is similar to such 
problems. For those problems, an attractor with a small 
attractor basin corresponds to a solution that is hard to 
find. Carrying the analogy further, we should not be sur- 
prised if a combinatorial problem that seem to be impos- 
sible to solve has plenty of solutions. Random Boolean 
networks could provide a playground for algorithms that 
handle such problems. 
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